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We present a new version of a dynamical spectral model for Large Eddy Simulation based on 
the Eddy Damped Quasi Normal Markovian approximation [l], 0] • Three distinct modifications are 
implemented and tested. On the one hand, whereas in current approaches, a Kolmogorov-like energy 
spectrum is usually assumed in order to evaluate the nonlocal transfer, in our method the energy 
spectrum of the subgrid scales adapts itself dynamically to the large-scale resolved spectrum; this 
first modification allows in particular for a better treatment of transient phases and instabilities, as 
shown on one specific example. Moreover, the model takes into account the phase relationships of 
the small-scales, embodied for example in strong localized structures such as vortex filaments. To 
that effect, phase information is implemented in the treatment of the so-called eddy noise in the 
closure model. Finally, we also consider the role that helical small scales may play in the evaluation 
of the transfer of energy and helicity, the two invariants of the primitive equations in the inviscid 
case; this leads as well to intrinsic variations in the development of helicity spectra. Therefore, our 
model allows for simulations of flows for a variety of circumstances and a priori at any given Reynolds 
number. Comparisons with Direct Numerical Simulations of the three-dimensional Navier-Stokes 
equation are performed on fluids driven by an ABC (Beltrami) flow which is a prototype of fully 
helical flows. Good agreements are obtained for physical and spectral behavior of the large scales. 



PACS numbers: 47.27.E-, 47.27.om, 47.27.ep, 47.27.er 



I. INTRODUCTION 

Turbulent flows are ubiquitous, and they are linked 
to many issues in the geosciences, as in meteorology, 
oceanography, climatology, ecology, solar-terrestrial in- 
teractions and fusion, as well as the generation and en- 
suing dynamics of magnetic fields in planets, stars and 
galaxies due to e.g. convcctive fluid motions. As manifes- 
tations of one of the last outstanding unsolved problems 
of classical physics, such flows form today the focus of 
numerous investigations. 

Natural flows arc often in a turbulent state driven by 
large scale forcing (novae explosions in the interstellar 
medium) or by instabilities (convection in the sun). Such 
flows involve a huge number of coupled modes at different 
scales leading to great complexity both in their temporal 
dynamics and in their emerging physical structures. Non- 
linearities prevail when the Reynolds number Rv - which 
measures the amount of active temporal or spatial scales 
in the problem - is large. In the Kolmogorov framework 
Q, the number of degrees of freedom increases as Rv 9//4 
for Rv ^ 1; for example in geophysical flows, Rv is often 
larger than 10 s . The ability to probe large Rv, and to 
examine in details the large-scale behavior of turbulent 
flows depends critically on the ability to resolve such a 
large number of spatial and temporal scales, or else to 
model them adequately. 

Only modest Reynolds numbers can be achieved by 
Direct Numerical Simulation (DNS) with nowadays com- 
puters. One way around this difficulty is to resort to 
Large Eddy Simulations (or LES, see e.g. H, 0] and 
references therein). Such techniques are widely used in 



engineering contexts, as well as in atmospheric sciences 
and, to a lesser extent, in geophysics (see Q) and astro- 
physics. Another class of models is based on two-point 
closures, like the Eddy Damped Quasi Normal Markovian 
approximation, or EDQNM [IJ. These models, developed 
in the mid seventies, gave rise to successful LES when 
taking their eddy viscosity formulation 0, . Such LES 
techniques have been used mostly in conjunction with 
pseudo-spectral methods, since being best expressed in 
Fourier space in terms of energy spectra. 

In this paper, we propose a new LES formulation that 
generalizes the usual EDQNM approach which is based 
on a Kolmogorov fc~ 5 / 3 spectrum (K41 hereafter), by al- 
lowing for a priori any kind of energy spectrum as may 
occur in the complex dynamical evolution of various tur- 
bulent flows, since our method (as will be shown later) is 
based on the evaluation of the transfer terms for energy 
and helicity. For example, there are small intermittency 
corrections to the K41 spectrum due to the presence of 
strong localized vortex filament structures in fluid tur- 
bulence; similarly, the presence of waves may alter the 
energy spectrum (see, e.g., OH). The method proposed 
here may also be particularly important when dealing 
with magnctohydrodynamics (MHD) flows, i.e. when 
coupling the velocity to the magnetic induction. In that 
case, the energy spectra can be either shallower 11| or 
steeper than fc -5 / 3 , because of anisotropy induced by a 
uniform magnetic field leading to Alfven wave propaga- 
tion and to weak turbulence for strong magnetic back- 
ground fl2| . Similarly, in the case of strong correlations 
between velocity and magnetic fields spectra that 
differ from the classical K41 phenomenology may emerge. 
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Note that, although we focus here on neutral flows, the 
extension of our model to conductive MHD flows presents 
no particular difficulties [T3 |. 



II. MODEL DESCRIPTION 



A. Equations 



In general, the traditional formulation of turbulent en- 
ergy transfers only takes into account the energy of the 
flow but not its helicity. However, the closure transfer 
terms for helicity are well known in the helical case @, 
including in MHD The kinetic helicity H = 1/2 < 
v ■ w > (where w = V x v is the vorticity) represents the 
lack of invariance of the flow by plane symmetry (w is an 
axial vector) . This global invariant of the Euler equation 
16 [has been little studied until recently (see however 
a 13, EH Eil and references therein). Furthermore, the 
intermittent structures that populate a turbulent flow at 
small scales, namely the vortex filaments, are known to 
be helical (see [HI, Ell, Ell); this implies that the nonlin- 
ear transfer terms involving small scales are weakened. 
This is consistent with several recent findings, namely: 
(i) helical vortex tubes, in a wavelet decomposition of a 
turbulent flow into a Gaussian component and a struc- 
ture component, represent close to 99% of the energy 
and corresponds to the strong tails of the probability 
distribution function of the velocity gradients [2(|; (ii) 
in a decomposition of the velocity field into large V and 
small v components, dropping (artificially) the nonlocal 
(in scale) nonlinear interactions (vV) leads to less inter- 
mittency [2lj . indicating that intermittency involves in- 
teractions between structures (like vortex tubes) that in- 
corporate small scales and large (integral) scales through 
a large aspect ratio; and (iii) the spectrum of helicity 
is close to fc -5 / 3 in the K41 range for energy, but not 
quite: the relative helicity 5 = H(k)/kE(k) decreases 
more slowly than 1 /k [22l . |23I | , indicating that the return 
to full isotropy is not as fast as one may have conjec- 
tured in the small scales. Finally, helicity is also invoked 
as possibly responsible for the so-called bottleneck effect, 
i.e. the accumulation of energy at the onset of the dis- 
sipation range [24[ , although it is not clear whether this 
effect is or not an inertial range phenomenon [23| . 



Let us consider the Fourier transform of the velocity 
v(x,i) and the vorticity w(x, t) = V x v(x, t) fields at 
wavevector k: 



/oo 
v(x, t) e - lk - x dx 
-OO 

/CO 
w(x,t)e- ikx dx. 
-OO 



(1) 



(2) 



In terms of the Fourier coefficients of the velocity compo- 
nents, the Navier-Stokes equation for an incompressible 
flow, with constant unit density, reads: 



(d t + isk 2 )v a (k,t) = £(k,t) + i^(M) 



(3) 



where F v is the driving force, v is the kinematic viscosity 
and t"(k, t) is a bilinear operator written as: 

tS(k,t) = -iP a/J (k)fc y Yl w/jkMKfa*) ; ( 4 ) 

p+q=k 

P a| a(k) = S a p — k a kp/k 2 is the projector on solenoidal 
vectors. In the absence of viscosity, both the total kinetic 
energy E = 1/2 < v 2 > and the total helicity H = 1/2 < 
v • w > are conserved and they are thus thought to play 
an important dynamical role in the temporal evolution 
of the fluid. The direct cascade of energy to the small 
scales and the related cascade of helicity [26| stem from 
these conservation laws. Furthermore, because helicity 
is thought to play an important role in the small scales 
(see e.g . | 16l| and references therein and more recently 
[la . [20 . [2J]), we are taking in this paper the approach 
of following the time evolution of both the energy and 
helicity spectra (see below). Taking the rotational of Eq. 
§5§ in Fourier space leads to: 



(d t + vk 2 )w a (k, t) = C(k, t) + F™(k, t) 



(5) 



Our dynamical spectral LES model, based on the 
EDQNM closure, is described in Section II. In Section 
III, numerical tests of the model are performed by com- 
parisons with three-dimensional direct numerical simu- 
lations (DNS) for strong helical ABC flows [25[, as well 
as a Chollet-Lesieur approach Predictions for high 
Reynolds number flows are also given. Section IV is the 
conclusion. Finally, details on closure expressions of the 
nonlinear transfers for energy and helicity, and on the 
numerical implementation of the model are respectively 
given in Appendix A and B. 



with 

COM) = e a sp k s k y ^ V 0(P> t ) v i( c L^ t ) > ( 6 ) 

p+q=k 

F™(k,t) = ie aS( 3 k s F£(k,i) . (7) 

We respectively define the modal spectra of energy 
£(k, t) and helicity 7i(k, t) in the usual way as: 



£(M) = iv(k,t).v*(k,t) - 



H(M) = -v(k,i)-w*(k,i) 



(8) 



(9) 
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where stars stand for complex conjugates. Note that w 
is a pseudo (axial) vector, and correspondingly the he- 
licity is a pseudo scalar. The integration of £(k, t) and 
H(k, t) over shells of radius k = |k| respectively gives 
the isotropic energy E(k,t) and helicity H(k,t) spec- 
tra. Their spatio-temporal evolutions obey the following 
equations: 

(d t + 2vk 2 )E{k 1 t) = T E (k,t) + F E (k,t) (10) 
{d t + 2vk 2 )H(k,t) = T H (k,t) + F H (k,t) (11) 

where T E (k,t) and Tji[k,t) denote energy and helicity 
nonlinear transfers at wave number k. They arc func- 
tional of the tensors involving triple correlations be- 
tween v(k,i), v(p,t) and v(q, t) with the constraint 
that p + q = k due to the convolution term in Fourier 
space emanating from the nonlinearities of the primitive 
Navier-Stokes equation. 

Finally, under the closure hypothesis customary to the 
EDQNM approach (see and references therein), the 
time evolution of E(k, t) and H(k, t) can be described by 
the EDQNM equations where the exact transfer terms 
T E (k,t) and Tn(k,t) in the equations above are re- 
placed by the closure evaluations denoted as T E (k, t) and 
Tn(k,t). The closure is done at the level of fourth-order 
correlators which are expressed in terms of third-order 
ones, with a proportionality coefficient - dimensionally, 
the inverse of a time - taken as the sum of all character- 
istic rates appearing in a given problem, namely the lin- 
ear (dispersive), nonlinear and dissipative rates. Tested 
against DNS Q , these closures allow for exponential dis- 
cretization in Fourier space and hence for exploration of 
high Reynolds number regimes. Their drawback is that 
all information above second-order moments is lost and 
phase information among Fourier modes is lost as well, 
so that, for example, intermittency is not present in this 
approach, nor are spatial structures. 

The full formulation of the EDQNM closure leads to a 
set of coupled intcgro-differcntial equations for the energy 
and helicity spectra E(k, t) and H(k, £), with the nonlin- 
ear transfers decomposed into emission terms (S El , S Es 
and , Sh 3 ), and absorption terms (S E2 , S Ei and Sh 2 > 
Shi)] note that we use S Ei and Shu with i £ [1,4], 
as short-hand notations for the full spectral functions 
S Ei (k,p, q,t) and (k,p, q, t). The expressions of these 
closure transfer terms are given in Appendix A. Note 
that absorption terms are linear in the spectra, the dy- 
namical evolution of which we are seeking, whereas emis- 
sion terms are inhomogeneous terms involving the p, q 
wavenumbers on which the double sum is taken (with 
p+q = k). The absorption term, S E2 , leads to the classi- 
cal concept of eddy viscosity, whereas the emission term, 
S El , is in general modeled as an eddy noise, although it 
is known through both experiments and DNS that the 
small scales are far from following a Gaussian distribu- 
tion, with substantial wings corresponding to strong lo- 
calized structures. Here, we present a different and novel 
method to treat the emission term. 



B. Spectral filtering 

When dealing with an LES method, as a complement 
in the unresolved small scales to the dynamical evolu- 
tion of the large scales following the Navier-Stokes equa- 
tion, we need to partition Fourier space into three re- 
gions. This means that we need to introduce a buffer 
region between the scales that are completely resolved 
(above k ~ , where k c is a cut-off wavenumber depend- 
ing on the resolution of the LES run), and the scales 
that are completely unresolved, say beyond ak c with a 
of 0(1). Following [III and according to the Test Field 
Model closure, the contribution of subgrid scales, to the 
cxplicitely resolved inertial scales, leads to an eddy vis- 
cosity depending both on the wavenumber and the en- 
ergy spectrum at that wavenumber. It is also shown that 
beyond 2k c , 85% of the transfer is covered by the eddy 
viscosity, while beyond 3k c , about 100% is covered; we 
thus choose to take a = 3. 

More specifically, the truncation of equations ([3]), ([5]), 
(fTQ)) and (jlip at two different wavenumbers k = k c and 
k = Sk c gives rise to three types of transfer terms, corre- 
sponding respectively to local, nonlocal and highly non- 
local interactions (where locality refers to Fourier space, 
i.e. interactions between modes of comparable wavenum- 
ber): 

(i) the fully resolved transfer terms Tg and Tjj involve 
triadic interactions such that k, p, and q are all three 
smaller than k c ; this interval is denoted A < ; 

(ii) the intermediate nonlocal tranfer terms Tg H , in 
which p and/or q are contained in the buffer zone be- 
tween k c and 3fc c (hereafter denoted A > ); and 

(iii) the highly nonlocal tranfer terms Tg#, in which p 
and/or q are larger than 3k c (hereafter denoted A>>). 

We choose to model Tg H and Tg-^ in Eqs. (JTUJ) and 
(fTTj) by appropriately modified EDQNM transfer terms. 
We therefore need to know the behavior of both energy 
and helicity spectra after the cut-off wavenumber k c = 
N/2 — 1, where N is the linear grid resolution of the 
numerical simulation. Whereas it is customary to assume 
a A: -5 / 3 Kolmogorov spectrum in this intermediate range, 
here we choose a different approach, namely, between 
k = k c and k = 3k c , both spectra are assumed to behave 
as power-laws (with unspecified spectral indices) followed 
by an exponential decrease, viz.: 

E(k,t) = E Q k- aE e- SEk , k c <k<3k c (12) 
H(k,t) = H k- aH e- SHk , k c < k < 3k c (13) 

where a E , S E , Eq, and ajj, Sh, Hq are evaluated, at 
each time step, by a mean square fit of the energy and 
helicity spectra, respectively. Note that it is understood 
that the Schwarz inequality |i2~(A;)| < kE(k) is fulfilled 
at all times. When either 8 E or Sh is close to zero, we 
consider that the energy (or helicity) spectrum has an 
infinite inertial range with a k~ aE H power law (see eq. 



4 



HH)), so we can write: 



E(k,t) 
H(k,t) 



H k~ aH , 



3k c < k < oo . 
3fc c < k < oo 



(14) 
(15) 



C. Eddy viscosity 

In the context of spectral models for the Navier-Stokes 
equation, the concept of eddy viscosity was introduced 
by Kraichnan [27| ■ This transport coefficient, denoted 
v(k\k c , i), allows to model the nonlinear transfer through 
a dissipative mechanism, as first hypothesized by Hciscn- 
berg. With and Tg > terms defined above, and where 
the hat denotes the fact that the EDQNM formulation of 
these partial transfers is taken, the eddy viscosity reads: 



v(k\k c ,t) 



T>(k,t)+T>>(k,t) 
2k 2 E{k,t) 
u > (k\k c ,t) + v»(k\k c ,t) , (16) 



thus separating the contribution stemming from the 
buffer zone (A>) and from the outer zone (A >> ). Note 
that only the part of the transfer proportional to the en- 
ergy spectrum at wavenumber k (i.e. the SE 2 {k,p,q,t) 
part defined in Appendix A) is taken into account in 
the derivation of v > {k\k Cl t). Indeed, in our model, the 
closure transfer term Tg(k,t) is integrated at each time 
step, but an eddy viscosity from this whole transfer term 
cannot be extracted; only the part that cxplicitcly con- 
tains E(k,t) (the linear part of the transfer in E(k,t)) 
enables the derivation of an eddy viscosity, namely: 



v>(k\k c ,t) 



° kpq SE 2 {k 1 p,q,t) 
A > 2k 2 E(k 1 t) 



dpdq 



A > kpq 2k 2 q 



[xy + z 3 )E(q,t)dpdq 



Let us now evaluate the eddy viscosity in the outer region, 
v >> (k\k c , t), coming from the highly nonlocal EDQNM 
transfer terms. Since the (k,p,q) triangles are very elon- 
gated in the A>> zone, with k <C p,q, an algebraical 
simplification occurs leading to an explicit expression for 
v >> (k\k c ,t). Indeed, it has been shown [2|| that a Tay- 
lor expansion of T^ > (k,t) with respect to k/q leads, at 
first order, to the following asymptotic transfer term: 



T^(k) 



-k 2 E(k) 
15 -'3k c 



S h „[5E(jp) 



8E(p) 
P ^p- ]dp 



(17) 

where time dependency is omitted for simplicity. Since 
now Tfi > (k) explicitcly depends on E(k). it is straight- 
forward to formulate the corresponding eddy viscosity 



'(k\k c ,t) thus defined as: 



v»{k\k et t) = 



15 



3Av 



dE(p) 



(18) 



When E(p) is replaced by its power law-exponential de- 
cay approximation (see Eq. I[14[) ). we recover the so- 
called "plateau-peak" model [!|: 



v»{k\k c ,t) 



1 — OiE 



a E C 



K 



E{3k c ,t) 



3kc 

(19) 

Finally, in the energy equation Eq. (TiT))) . the total eddy 
viscosity derived from the nonlocal and highly nonlocal 
transfer terms is simply obtained by adding the two con- 
tributions, as stated before: v(k\k c ,t) = v > (k\k Cl t) + 
v»(k\k c ,t). 

Note that, in the helicity equation Eq. (jXTJ) , the trans- 
port coefficient stemming from the helicity transfer term 
Th (k, t) can be similarly evaluated in the buffer zone and 
the outer zone, and written as: 



u H (k\k c ,t) 



T>{k,t)+T>>{k,t) 
2k 2 H(k,t) 
iy>(k\k c ,t) + iy>>(k\k c ,t), (20) 



with 



u>(k\k c ,t) 



kp SH 2 (k,p,q,t) 



dpdq. (21) 



, A > 2k 2 H(k,t) 

It is straighfoward to show that this eddy viscosity part 
has the same formulation than v > (k\k C: t) (see Appendix 
A for SH 2 (k,p,q,t) definition). For the helicity transfer 
term T^ > (k,t), simple algebraic calculations lead to: 



T. 



h (k) = --k 2 H(k) 



3k c 



9 kpp [5E(p)+p^^}dp. 



dp 

(22) 

The integrands in Eqs. (|17|) and (f2"2")) are thus identi- 
cal; this in turn provides the same eddy viscosity in the 
outer domain than for the energy, namely {k\k c , t) = 
v >> (k\k c ,t). Altogether, the same total eddy viscosity 
appears in both the energy and helicity equations. This 
is expected from the formulation of the spectral closure, 
in which the temporal dynamics of the second-order ve- 
locity correlation function is separated into its symmetric 
(energetic) and anti-symmetric (helical) parts. 



D. Helical eddy diffusivity 

At wavenumber fc, the energy transfer obtained from 
the use of the EDQNM closure involves a linear term in 
the helicity spectrum H(k,t) (specifically, SE 4 (k,p,q,t) 
defined in Appendix A, Eq. IA5[) : from this term, a new 
transport coefficient, similar to the v'' (k\k c ,t) eddy vis- 
cosity, can be built. In the buffer zone, this new coeffi- 
cient, hereafter named "helical eddy diffusivity", reads: 

8 kpq SE 4 (k,P, q, t) 



v>{k\k c ,t) 



A> 



'kpq 



2k 2 H(k,t) - dpdq 

z(l-y 2 )H(q,t)dpdq 



2k 2 q 



(23) 
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Note that, dimcnsionally, this helical diffusivity v scales 
as vjk. As before a total helical eddy diffusivity can be 



eddy diffusivity, namely: 



defined as v{k\k c ,t) = v > (k\k c ,t) + v >y (k\k c ,t), where (dt + vk 2 )v a (k, t) = — iPo,^(k)fc 7 vp(p,t)v 1 (q) 



v{k\k Cl t) represents the contribution of the small-scale 
helicity spectrum to the kinetic energy dissipation. Re- 
call that, in the outer zone, the Taylor expansion of 
the highly nonlocal transfer, (k,t), with respect to 
k/q <C 1, leads at first order to Eq. (JT7J) , with no lin- 
ear contribution from H(k,t). We therefore assume that 
the transfer part associated with helical motions in the 
outer zone is negligible, such as v >> (k\k c ,t) = 0. The 
total helical eddy diffusivity thus reduces to V(k\k c , t) = 
V>(k\k c ,t). 



E. Emission transfer terms 

The parts of the EDQNM transfer terms which arc 
not included either in the eddy viscosity or in the helical 
eddy diffusivity, involve energy and helicity interactions 
at wavenumbcrs p and q both larger than k c . Respec- 
tively denoted T/ 9 (fc,t) and f^ q (k,t), they read: 



f/*(M) 
fp(k,t) 



3fc c pk+p 



k c J k—p 
3k c pk+p 



kpa (*)(&i + S E3 )dpdq 



k c J k—p 



°k Pq ( t ){S Hl + S H3 )dpdq 



(24) 



where SE it Hi stands for SEi.Hi(k,p,q,t). 

On the one hand, the established eddy viscosity and 
helical eddy diffusivity can be directly used in the Navier- 
Stokes equation for the modal energy and helicity spec- 
tra, £(k,t) and H(k,t) respectively. On the other hand, 
in order to implement in these modal equations, the 
isotropic transfers T^ q (k,t) and T^ q (k,t), we assume 
that they are uniformly distributed among all k wavevec- 
tors belonging to the same k-shell. This means that the 
nonlocal modal energy and helicity transfers, respectively 
T^ q (k,t) and T£ q (k,t), can be expressed as T^ q (k,t) = 
T p E q (k,t)/Aitk 2 andfg 9 (M) =fP q (k,t)/4irk 2 . 



F. Numerical field reconstruction 

To compute our LES model for all k < k c , we proceed 
in two steps. At a given time, the Navier-Stokcs equation 
is first solved using the eddy viscosity and the helical 



p+q=k 

k,p,q<k c 

—v(k\k C) t)k 2 v a (k, t) 
— v(k\k c , t)k 2 w a (k, t) 
+K(Kt). (25) 

Then, the effects of the emission terms, T E pq (k,t) and 
Tjf q (k,t), are introduced in the numerical scheme. In 
most previous studies, these terms are taken into ac- 
count through a random force, uncorrelated in time (see 
e.g. HH); this corresponds to the vision that they rep- 
resent an eddy noise originating from the small scales. 
However, the small scales are all but uncorrelated noise; 
the phase relationships within the small-scale structures 
play an important role, albeit not fully understood, in 
the flow dynamics. It is well-known that a random field 
with a fc -5 / 3 Kolomogorov energy spectrum, but other- 
wise random phases of the Fourier coefficients, is very 
different from an actual turbulent flow, lacking, in par- 
ticular, the strong vortex tubes so prevalent in highly 
turbulent flows. Similarly, it has been recently shown 
[30j | that, upon phase randomization, the ratio of non- 
local energy transfer (i.e. the transfer involving widely 
separated scales) to total energy transfer reduces to a 
negligible amount, whereas this ratio is close to 20 % at 
the resolutions of the performed numerical experiments, 
corresponding to a Taylor Reynolds number of about 10 3 . 
These considerations lead us to directly incorporate the 
emission terms in the second step of the our numerical 
procedure. The modal spectra of the energy and the he- 
licity, associated to the v(k,t) field computed from Eq. 
([2~5"T) , now has to verify the following equations where the 
emission transfer terms are taken into account; 



(d t + 2vk 2 )£{k,t) = 



(d t + 2vk 2 )H{k, t) 



-2v(k\k c ,t)k 2 £(k, t) 
-2v(k\k c ,t)k 2 H(k, t) 

+^<<M) + %2 

+F E {k, t), 

-2v(k\k c ,t)k 2 H(k,t) 
-2v(k\k c ,t)k 4 £{k, t) 
T P H g (k,t) 



(26) 



+T<(k,i) 
+F H (k,t) 



47rfc 2 



(27) 



where Te (k, t) and Th (k, i) denote the spectral terms 
stemming from the driving force. Recall that T^(k,t) 
and Tfl(k, t) are the resolved transfer terms based on 
triadic velocity interactions with k, p, and q all smaller 
than k c . Once the uptaded £(k, t) and H(k 7 t) modal 
spectra are obtained, the velocity field is updated. How- 
ever, a difficulty immediately arises: the phase rela- 
tionships between the three components of the velocity 
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field in the EDQNM (and other) closures is of course 
a priori lost. We thus proceed to the reconstruction 
of the three spectral velocity components, written as 
v a (k,t) = p a (k, t)e ll ^ Q ( k '*), and rebuild the different ve- 
locity phases by using the incompressibility and realis- 
ability (|W(k, i)| < k£(k, t)) conditions, as explained in 
Appendix B. 



III. NUMERICAL TESTS OF THE MODEL 



A. Numerical setup 

In order to assess the model accuracy to reproduce the 
physics involved in fluid flows, we performed Direct Nu- 
merical Simulations (DNS) of the Navier-Stokes equation 
and computations using our LES formulation. We de- 
note LES P the code with partial recovery of phases and 
without helical effects (i.e., with v = and = 0), and 
LES PH the code with helical effects incorporated. In our 
LES description, the energy spectra - and helicity spectra 
when considered - of the subgrid scales self adapt to the 
large scale resolved spectra {i.e. no spectral scaling laws 
are prescribed). We can therefore study a variety of flows, 
such as either low or high Reynolds number flows, or the 
early phases of the temporal development of flows when 
Kolmogorov spectra are not yet established. Indeed, the 
cut-off wavenumber, k c , can as well lie in the dissipa- 
tion range instead of the inertial range which is the case 
of standard LES approaches based on closures together 
with a fc -5 / 3 Kolmogorov spectrum. This enables accu- 
rate comparisons with DNS at a given viscosity. We also 
compared our numerical approach to a Chollct-Lcsicur 
LES model (CL) 0, where a v = 2.e kinematic vis- 
cosity is added to the turbulent viscosity for comparison 
purpose (see Table [J. The codes use a pseudo-spectral 
Fourier method in a [0— 27r] 3 -periodic box and an Adams- 
Bashforth second-order scheme in time. 

To test the ability of our LES models to reproduce heli- 
cal flow features, we focus on flows driven by a prototype 
Beltrami flow (v = ±w), namely the ABC flow (see e.g. 
H): 



FABc(fco) 



B cos(fco2/) + C sin(fco^) 
C cos(k z) + A sm(k x) 
Acos(k x) + Bsm(k n y) 



(28) 



with k = 2, and A = B = C = 1. 

The run parameters are summarized in Table [U The 
definitions used for the integral scale, L, and the Taylor 
microscale, A, are based on the kinetic energy spectrum 
E(k); respectively L = 2tt Jk~ 1 E(k)dk/ J E(k)dk and 

A = 2tt [/ E(k)dk/ J k 2 E(k)dk] 1/2 . Note that the char- 
acteristic flow quantities - L and A scales, r.m.s. velocity, 
U rms , nonlinear time scale, t nl , and Reynolds number, 
Rv - are time averaged quantities once the steady state is 
achieved for the computed flow. 



B. Spectral features 

We first investigate the flow spectral behavior on one- 
dimensional energy, enstrophy and helicity spectra ob- 
tained from the different models. These spectra are aver- 
aged over 67 nonlinear turnover times spanning the flow 
steady phase from t = 10.0 up to t = 60.0, the final 
time reached in the simulation. Fig. [T] shows the time 




FIG. 1: Time averaged energy spectra < E(k, i) > for data I 
(256 3 DNS, solid line), II (64 3 LES PH, dashed line) and III 
(64 3 LES P, dotted line). See Table H 

averaged energy spectra < E(k,t) > for LES and DNS 
computations (runs I, II and III). Both LES results show 
good agreement with the corresponding DNS ones, up to 
k c — 31, the maximum wavenumber of the LES calcu- 
lations. Small differences are observed in both LES P 
and LES PH models at the largest wavenumbers. When 
looking at the vorticity density spectra, these differences 
are amplified as small scales arc emphasized (see Fig. 
However, the mean characteristic wavenumber of the ve- 
locity gradients, defined as the maximum of the time av- 
eraged vorticity density spectrum, corresponds to k = 9 
in all runs. 

Note that, for the different flows, the time averaged 
Reynolds number and characteristics integral scale are al- 
most the same (see Table HI, while the Taylor microscale, 
A, and its associated Reynolds number, R\, differ due 
to non negligible intensities of the enstrophy spectrum 
(namely, k 2 E(k) used to compute A, with isotropy as- 
sumed) after the cut-off wavenumber k c . When down- 
sizing the 256 3 DNS data to 64 3 grid points by filtering 
all wave vectors k such as |k| > k c (case Ir in Table HJ, 
the Taylor small scale quantities obtained from Ir, LES 
PH II and LES P III runs get closer. Thus, our LES 
models can estimate the mean correlation length scale of 
the vorticity when the smallest scales are properly filtered 
out. 

The time averaged helicity spectra are plotted in Fig. [3] 
At large scale, the LES P and LES PH models provide 
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TABLE I: Parameters of the simulations I to IX. Linear grid resolution N, kinematic viscosity v and time averaged quantities: 
Taylor microscale A and integral scale L\ r.m.s.velocity U rma =< v 2 > 1,/2 ; integral Reynolds number Rv = U rms L/v; eddy 
turnover time tnl = L/U r ms ; Im is the final time of integration. Note that the Ir label stands for 64 3 reduced data obtained from 
the 256 3 DNS computation. The LES P (vs. PH) label stands for our model computations without (resp. with) incorporating 
the helicity transport coefficients v and emission transfer terms TjJ- 9 . The LES CL label stands for a Chollet-Lesieur scheme 
where the kinematic viscosity is added to the scheme eddy viscosity (see text). 







N 


V 


A 


L 


U rms 


Rv 


T N L 


£m 


I 


DNS 


256 


5.e~ 3 


0.81 


2.38 


3.19 


1525 


0.75 


60 


Ir 


Reduced DNS 


64 


5.e" 3 


0.92 


2.38 


3.19 


1530 


0.75 


60 


II 


LES PH 


64 


5.e~ 3 


0.93 


2.37 


3.20 


1519 


0.74 


60 


III 


LES P 


64 


5.e" 3 


0.93 


2.40 


3.22 


1544 


0.74 


60 


IV 


DNS 


512 


2.e" 3 


0.49 


2.32 


3.34 


3881 


0.70 


7.0 


V 


LES PH 


128 


2.e" 3 


0.59 


2.30 


3.36 


3877 


0.69 


7.0 


VI 


LES P 


128 


2.e~ 3 


0.59 


2.33 


3.37 


3925 


0.69 


7.0 


VII 


LES CL 


128 


2.e" 3 


0.66 


2.38 


3.29 


3928 


0.72 


7.0 


VIII 


LES PH 


256 


5.e- 4 


0.36 


2.47 


3.38 


16693 


0.73 


10 


IX 


LES P 


256 


5.e" 4 


0.36 


2.35 


3.31 


15565 


0.71 


10 




k 

FIG. 2: Time averaged enstrophy spectra for runs I (solid 
line), II (dashed line) and III (dotted line). 
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FIG. 3: Time averaged helicity spectra < H(k,t) > for runs 
I (solid line), II (dashed line) and III (dotted line). 



a close approximation of the DNS helicity spectrum up 
to k ~ 20. For k > 20, the LES PH model, designed 
to take into account helical effects, slightly overestimates 
the < H(k, t) > magnitudes of DNS data, while the LES 
P model dissipates too much helicity. 



C. Temporal evolution 

In this section, to study the temporal behavior of the 
flows, in the spirit of the analysis performed for freely 
evolving fluids, we focus on the temporal phase before 
the steady state regime. Fig. 0] shows the evolution of 
the global kinetic energy E(t) and helicity H(t) for DNS 
together with LES computations, namely runs IV (512 3 
DNS), V (128 3 LES PH) and VI (128 3 LES P). We first 
observe that for both LES models, with and without heli- 



cal effects, energies closely follow the growth phase of the 
DNS energy. Indeed, during the inviscid phase, t < 1.0, 
the small scales are generated with negligeblc effects on 
large scales, since their intensities are very weak. Thus, 
at that times, the LES modelling has only a reduced ac- 
tion. Later on, during the following growth phase, up to 
t ~ 2.3, the effect of the subgrid scales onto the resolved 
ones becomes important, and the LES models correctly 
reproduce the DNS dynamics again. Differences then 
start to appear between the LES energy approximations 
and the DNS energy, as all scale intensities increase, and 
therefore so does the influence of the intermediate and 
highly nonlocal transfer terms (see Section II. B). How- 
ever, their mean values stay close to the DNS energy 
(see U rms in Table [TJ) . The same remarks hold for the 
temporal evolution of the kinetic helicity. However, in 
this latter phase, t > 5, one can note that the LES PH 
model provides a slightly better approximation than the 
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The temporal behavior of the LES flows can be un- 
derstood when looking at kinetic energy spectra at early 
times, as plotted in Fig. El Instantaneous DNS energy 
spectra are well fitted by LES spectra up to t ~ 3.0, 
including in the phase of development toward a Kol- 
mogorov spectrum. At larger times, the instantaneous 
modeled slightly spectra diverge from the DNS ones 
, as seen in the steady state in the previous section, 
with small scales slightly underestimated and large scales 
slightly overestimated. 



FIG. 4: Evolution of kinetic energy E(t), lower curves, and he- 
licity H(t), upper curves, for runs IV (512 3 DNS, solid line), 
V (128 3 LES PH, plusses) and VI (128 3 LES P, triangles). 

LES P one, for both energy and hclicity. When comput- 
ing the temporal mean of the relative error between the 
LES and DNS data, we obtain for the energy 1.28% for 
the helical model (versus 1.36% for LES P). For the hc- 
licity, these errors are respectively 2.20% for LES PH and 
2.27% for LES P. Considering that the cost of computing 
the additional helical term is rather small (the LES PH 
needs 6% CPU time more than the LES P simulation), 
the slight improvement when using the helical model is 
worth considering; in particular, note that it reproduces 
better the temporal oscillatory variation of the total en- 
ergy, although at a higher intensity. On the other hand, 
the fact that the non-helical model performs almost as 
well shows that helicity does not play a significant role in 
the small scales, in agreement with the statistical argu- 
ment of return to isotropy in the small scales, and with 
the fact that the relative helicity decays faster than fc _1 
in the small scales which are being truncated in an LES 
computation. 




k 



FIG. 5: Temporal development of energy spectra E(k, t) 
shown at t — 0.5, t — 1.0, t = 2.0 and t = 3.0, for runs 
IV (DNS, solid line) and V (LES-PH, plusses). 



D. Statistical analysis 

Flow statistics are now investigated. Probability den- 
sity functions (hereafter, pdf), are computed from 20 ve- 
locity snapshots extracted each 3tjvl and spanning the 
flow steady states, between t = 10.0 and t = 60.0. We 
compare data sets obtained from runs I (256 3 DNS), Ir 
(64 3 reduced DNS data) and II (64 3 LES PH), for which 
we have large velocity samples. For clarity purpose, since 
differences between helical and non helical models are not 
visible on the pdf, we only represent LES PH results ver- 
sus DNS ones. Fig. [5] (a) displays the statistical distribu- 
tion of the Vz-velocity component, after being normalized 
so that a 1 —< >= 1 (a being the standard deviation), 
together with a Gaussian distribution. The obtained dis- 
tributions for all data runs are close to Gaussian, typical 
of large scale velocity behavior. The LES models being 
designed for recovering correctly the large-scale flow, the 
LES pdf are identical to those of the truncated DNS data 
set Ir (where the smallest DNS scales have been filtered 
out for |k| > k c ). Note that they are also very close to 
the velocity distribution of the full DNS data set. 

Examples of spatial distributions of longitudinal and 
lateral velocity derivatives, dv x /dx and dv y /dx respec- 
tively, are shown in Fig. [6](b) and Fig.[6](c). The distri- 
butions are closer to an exponential than to a Gaussian. 
This behavior is even more pronounced for lateral deriva- 
tives with a slight departure from an exponential law. 
The wings of the pdf are mainly due to small-scale veloc- 
ity gradients. Since the DNS flow has more excited small 
scales, the wings of the associated velocity derivatives 
distributions are more extended than for LES data. The 
LES distributions correctly reproduced the DNS ones up 
to 3cr, however there is almost no differences with the pdf 
of the Ir data set. 

TABLE II: Temporal mean skewness of velocity derivatives for 
runs I-III and Ir data (see Table U}. Error bars are computed 
from instantaneous data. 

I DNS Ir II LES PH III LES P 

dv x /dx -0.45 ± 0.05 -0.35 ± 0.05 -0.35 ± 0.04 -0.33 ± 0.05 
dvy/dy -0.45 ±0.06 -0.34 ± 0.04 -0.34 ± 0.06 -0.34 ± 0.06 
dv z /dz -0.46 ±0.06 -0.35 ± 0.05 -0.34 ± 0.07 -0.33 ± 0.05 
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FIG. 6: Mean probability distributions of v x (a), dv x /dx (b) 
and dvy/dx (c), normalized so that a = 1, for data I (256 3 
DNS, in blue), Ir (64 3 reduced DNS data, in black) and II 
(64 3 LES PH, in red), shown together with a Gaussian distri- 
bution (dotted line). See Table U 



TABLE III: Temporal mean flatness of the velocity gradients 
for the same runs as in Table [TT1 





I DNS 


Ir 


II LES PH 


III LES P 


dv x /dx 


5.0 ±0.2 


4.0 ±0.2 


4.0 ±0.3 


4.0 ±0.3 


dv y jdx 


7.2 ±0.6 


5.0 ±0.4 


4.9 ±0.3 


4.8 ±0.4 


dvz/dx 


7.3 ±0.4 


5.1 ±0.4 


4.9 ±0.6 


4.8 ±0.5 


dv x /dy 


7.4 ±0.4 


5.0 ±0.3 


4.9 ±0.5 


4.8 ±0.5 


dvy/dy 


5.1 ±0.2 


3.9 ±0.1 


3.9 ±0.2 


3.9 ±0.2 


dv z /dy 


7.3 ±0.6 


5.0 ±0.3 


4.8 ±0.4 


4.7 ±0.5 


dv x /dz 


7.3 ±0.6 


5.0 ±0.5 


4.9 ±0.5 


4.7 ±0.4 


dvy/dz 


7.3 ±0.9 


5.0 ±0.4 


4.8 ±0.7 


4.8 ±0.8 


dv z /dz 


5.1 ±0.4 


4.0 ±0.2 


3.9 ±0.3 


3.9 ±0.3 



In order to quantify the distributions of the velocity 
fluctuations, and their differences among DNS, LES PH 
and LES P data, we compute low-order moments, namely 
the skewness (S3) and flatness (S4) factors of the velocity 
derivatives, defined as S„ =< /" > / < f 2 >™/ 2 where 
/ stands for any velocity derivatives. Their temporal 
means and error bars, are given in Table ITT1 and Table iLUl 
respectively. For a fair comparison with the results of 
our LES simulations, we also compute the skweness and 
flatness factors based on the Ir reduced DNS velocity 
fields. The pdf of the longitudinal velocity derivatives 
present an asymmetry that yields their well-known neg- 
ative skweness, with S3 ~ —0.45 for DNS velocities, a 
value comparable to other simulations at R\ ~ 500 (e.g. 
HH). The skweness for the reduced DNS data (from 256 3 
to 64 3 ) is about 22% lower, a reduction simply due to the 
truncation of velocity fields in Fourier space as noted in 
[32^ . The LES models give almost identical results, with 
in most cases, slightly closer values for LES PH data sets. 
It remains an open problem to know whether this type 
of agreement persists for higher Reynolds numbers. For 
all runs, the lateral velocity derivatives arc much more 
symmetric: S3 ~ (with fluctuations varying from 10~ 2 
to 10~ 4 ), as expected for fields that are almost statis- 
cally isotropic [33| • The longitudinal and lateral velocity 
derivatives do not have the same S4 flatness factors. For 
the DNS data, the flatness values for dv x /dx, dv y /dy and 
dv z /dz are close to 5, while the flatness for the lateral 
velocity derivatives (for example dv x /dy and dv x /dz) 
are larger, with values around 7. The reduced Ir data 
present a loss of ~ 20% and ~ 30% for longitudinal and 
lateral derivatives, respectively. Once again, both LES 
data provide similar values with a better approximation 
of the lateral derivatives for the LES PH computation. 



E. Visualization in physical space 

The topological properties of the different flows are 
now investigated: comparisons between DNS and LES 
PH computations are carried out on either instantaneous 
and or mean velocity fields. 




FIG. 7: Contour plots of the velocity intensity | v(x,t) | at 
time t = 2.8 on the boundaries of the periodic box, for run 
IV (512 3 DNS, left) and run V (128 3 LES PH, right). 
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For runs IV (512 3 DNS) and V (128 3 LES PH), Fig. 
displays contour plots of the velocity intensity at time 
t = 2.8, shown on three sides of the periodic box. Al- 
though our LES model cannot exactly reproduce the 
DNS flow, one can notice that the main flow structures, 
and their intensities, are well represented at times before 
the statisticaly stationary regime. More precisely, the 
mean spatial correlation of the pointwise LES PH veloc- 
ity field with the DNS one is 84.37%, while it is 84.32% 
in the case with the LES P pointwise velocity at the same 
time (not shown). 




at t = 2.8 for data IV (512 3 DNS) and V (128 3 LES PH) 
(Fig. [7]), while for the mean flow shown in (Fig. [HJ, the 
coefficients arc based on the time averaged fields of runs 
I (256 3 DNS) and II (64 3 LES PH). For comparison, the 
isotropy coefficients are also computed with the data of 
the LES P runs VI (256 3 ) and III (64 3 ). One can see 
that, altogether, the isotropic properties of the flow are 
correctly restaured by the LES models. 

TABLE IV: Isotropy coefficients of velocity, CJ so , and vortic- 
ity fields, C™ D , for flows shown in (Fig. fTJ and (Fig. IS}. For 
runs I to III, values are computed from the mean flows. For 
runs IV to VI, they are given at t = 2.8. For completeness, 
the coefficient values are also given for LES P runs. 







C v 

^ ISO 


^iso 


I 


DNS 


0.990 


1.010 


II 


LES PH 


1.026 


0.979 


III 


LES P 


0.988 


1.009 


IV 


DNS 


0.961 


1.033 


V 


LES PH 


0.955 


1.035 


VI 


LES P 


0.955 


1.029 



FIG. 8: Isosurface of the mean flow intensity |< v(x,i) >| for 
run I (256 3 DNS, left) and run II (64 3 LES PH, right) plotted 
at a level of 2, with maximum intensity values of 3.18 for run I 
and 3.31 for run II. Averages are taken over 50 computational 
times spanning the steady states. 

For runs I (256 3 DNS) and II (64 3 LES PH), Fig. M 
shows one isosurface of mean velocity intensities at a level 
of 2, i.e. roughly at 63% of the maximum intensity for 
DNS data, and at 60% for LES PH data. The mean 
velocity fields are time averaged each 3tnl during the 
steady states that correspond to 67tn l for the DNS flow 
and 68ttvl for the LES PH run. From these runs, with 
the longest stationary phases, we observe that the mean 
velocity field, and thus the large and intermediate flow 
scales, are not alterated by the modeling of the transfers 
linked to the subgrid scales. Indeed, the mean spatial 
correlation between the pointwise LES and DNS time 
averaged velocity fields is 91.96% (for completness, it is 
92.38% when the LES P mean field is considered). 

Finally, the flow isotropy is estimated by means of 
coefficients computed as in [34| . For each wavevector 
k, an orthonormal reference frame is defined as (k/|k|, 
ei(k)/|ei(k)|, e 2 (k)/|e 2 (k)|), with ei (k) = k x z and 
e 2 (k) =kx ei(k), where z is the vertical unit wavevec- 
tor. In that frame, since the incompressibility condi- 
tion, say for the velocity field, yields k ■ v(k) = 0, 
v(k) is only determined by its two components Vi(k) 
and v 2 (k). The isotropy coefficient is then defined as 
CJ so = < |vi | 2 > / < |v 2 | 2 >, with thus a unit value 
for fully isotropic flows. A similar coefficient can be based 
on the vorticity field, C™f D , characterizing the small scale 
isotropy. Isotropy coefficient values are given in Tabic ITVl 
for velocity and vorticity fields of the flows visualized in 
(Fig. [7]) and (Fig. [5]) . Instantaneous values are computed 



F. Comparison with a Chollet-Lesieur approach 

It may be instructive to test our model against an- 
other spectral LES approach, also based on the EDQNM 
closure; we have thus performed a simulation using a 
Chollet-Lesieur scheme Q (run VII LES CL in Table H}. 
The CL model allows energy tranfer from subgrid to re- 
solved scales through a dissipation mechanism, with the 
help of a dynamical eddy viscosity vchik-,^) defined as : 



v CL {k,t) = Cv+{k,t)^E(k cut7 t)/k c 



(29) 



k C ut = N/2 — 3 is the cut-off wavenumber, N being 
the number of grid points per direction, and v + (k,t) 
is the so-called cusp function evaluated as v + (k,t) = 
(1 + 3.58(fc/fc cut ) 8 ). We recall that Cy/E(k cut ,t)/k cu t is 
the asymptotic expression of the nonlocal tranfers from 
subgrid to resolved scales, and it assumes a A: -5 / 3 Kol- 
mogorov spectrum extending to infinity. The constant 
C is adjusted with the Kolmogorov constant computed 
from the ABC flow resolved by the DNS run using 512 3 
grid points, namely C = 0.14. To be able to compare our 
DNS and LES simulations, runs IV to VI (see Table 0}, 
to a CL simulation with 128 3 grid points, the kinematic 
viscosity used for the former runs, v — 2.e -3 , is added 
to i>cL{k,t). The asymptotic value of the eddy viscos- 
ity can be obtained as the temporal mean of ^cl(Oj^) 
in the time interval [3.5, 7.0] and is here estimated to be 
~6.e- 4 . 

Both CL and LES PH runs provide a close agreement 
with DNS temporal evolutions of kinetic energy and he- 
licity during the growth phase, i.e up to t ~ 2.3 (see 
Fig. 0. In steady states, although noticeable deviations 
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FIG. 9: Evolution of energy E(t), lower curves, and helicity 
H(t), upper curves, for data IV (512 3 DNS, solid line), V 
(128 3 LES PH, plusses) and VII (128 3 LES CL, squares). 



occur between DNS and LES approximations, a slightly 
better assessment is visible for the LES PH model; in- 
deed, the peak of oscillation at t ~ 5.4 is weaker in run 
VII. Note that comparisons of LES PH versus LES P 
data, runs V and VI respectively, are already presented 
in Section III.C. 
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FIG. 10: Mean energy (top) and helicity (bottom) spectra 
for runs IV (512 3 DNS, solid line), V (128 3 LES PH, dashed 
line) and VII (128 3 LES CL, dotted line). 

Mean energy and helicity spectra, time averaged from 
t = 3.5 to t = 7.0, are plotted in Fig. QJ5] for runs IV, 
V and VII described in Table [J Energy spectra only 
exhibit a short k~ 5 ' 3 inertial range for all runs, including 
the CL run, due to the procedure we chose avoiding a 
non-zero asymptotic viscosity. With v = in the CL 
scheme, the results might be different. However, for both 
energy and helicity spectra, LES PH data give slightly 



closer results when compared to DNS data, as our CL 
simulation seems to overestimate positive energy transfer 
from resolved scales (between k ~ 15 to k cut = 61) to 
subgrid scales. 



G. Predictions for high Reynolds number flow 

In this last section, we present model computations for 
flows at high Reynolds number. Recently, Kuricn et al. 
[24j showed that for flows with maximum helicity, both 
energy an helicity spectra exhibit a fc~ 4 / 3 scaling range 
following the /s -5 / 3 Kolmogorov range and preceding the 
dissipation range, a result also found in [23|. This change 
in the energy spectrum is estimated from energy flux 
based on a characteristic time scale, denoted th, of dis- 
tortion (or shear) of eddies with wavenumber k submitted 
to out-of-plane velocity correlations, corresponding to he- 
licity transfer. We recall that the two dynamical times 
in competition are estimated by T H (k) ~ (\H(k)\k 2 /2) _1 
and T% L (k) ~ (£(fc)fc 3 ) -1 . Moreover, from DNS of the 
forced Navier-Stokes equation, these authors associate 
the well-known bottleneck effect with this scaling change 
when th becomes physically relevant. The ABC flow be- 
ing known for the presence of strong helical structures, we 
performed two simulations using our LES PH and LES P 
models at kinematic viscosities v = 5.e -4 , with 256 3 grid 
points (respectively run VIII and run IX in TableQJ . For 
these flows, the total helicity H(t) = 1/2 < v(t) -w(t) >, 
averaged over ~ 8tnl in the steady state, is equal to 
9.26 for data set VIII and to 8.65 for run IX. Tempo- 
ral means of the total relative helicity H(t)/E{t) are also 
close (within 1%) for both computations, namely 0.195 
for the LES PH run versus 0.186 for the LES P one. More 
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FIG. 11: Time averaged relative helicity < \H(k)\/kE(k) > 
for data sets : XIII (256 3 LES PH), solid line, and IX (256 3 
LES P), dashed line. A fc _1 slope is plotted for comparison 
(dotted line). 

precisely, in the range 10 < k < 100, the relative helic- 
ity \H(k)\/kE(k), viewed as an estimation of the ratio 
T H (k)/T^f L (k), lies in between 13.5% and 3.5% for run 
VIII, and falls from about 14% to 2% for run IX (see 
Fig. ITT]) , a typical ratio for strong helical flows [24| • Note 
also that the relative helicity obtained by both model 
scales closely to a k~ x power law observed in previous 
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DNS experiments [23[. 
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FIG. 12: Temporal mean of energy spectra compensated by 
k 5/3 (solid lines) and fc 4/3 (dash lines) for runs XIII (256 3 
LES PH) and IX (256 3 LES P) (see insert). Horizontal seg- 
ments indicate ranges of k~ 5 ^ 3 and k~ A ^ 3 scaling regime. 

Fig. [12] displays mean energy spectra for the two mod- 
eled flows compensated by fc 5//3 and fc 4 / 3 respectively. A 
fc -5 / 3 scaling appears in the range 4 < k < 10, followed 
by a fc~ 4 / 3 regime in for 10 < k < 40, and with no ap- 
pearance of a bottleneck effect, although, in the latter 
wavenumber interval the estimated time ratio th/t^j, 
ranges from 37% to 20%. Similarly, time averaged he- 
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FIG. 13: Temporal mean helicity spectra compensated by 
k 5 ^ 3 and k A ^ 3 for the same runs as in Fig. 1121 

licity spectra compensated by fc 5 / 3 and fc 4 / 3 are plotted 
in Fig. [13] A fc -5 / 3 behavior is seen in approximatively 
the same range than for energy spectra, while the £;~ 4 / 3 
scaling occurs from k ~ 10 to k ~ 60 for the LES P flow 
and up to k c , the maximun computational wavenumber, 
for the LES PH flow. Recall that, at high wavenumbcrs. 
our LES PH model slightly overestimates helicity spectra 
while our LES P model underestimates them. However, 
both LES models well reproduce the observed spectral 
behaviors obtained from DNS of flows at lower kinematic 
viscosities iv = l.e~ 4 and v = 0.35e~ 4 ) in [24| . 



IV. CONCLUSION 

In this paper, we derive a consistent numerical method, 
based on the EDQNM closure, to model energy inter- 



actions between large and small scales for the Navier- 
Stokes equation. As no spectral behavior is a priori 
given, our dynamical LES method allows for the mod- 
cling various flows, whether turbulent or not, compared 
to former spectral models which can, in principle, sim- 
ulate only infinite Reynolds number flows. The phase 
relationships of the small scales are taken into account 
through a numerical reconstruction scheme for the spec- 
tral velocity field. Helical effects in turbulcnts flows, like 
in vortex filaments, are also considered through the eval- 
uation of the energy and helicity transfers. For this pur- 
pose, an " helical eddy diffusivity" , similar to an eddy 
viscosity, and the emission transfer terms in which the 
helicity spectrum appears, are incorporated in a second 
model. Numerical tests of our two methods, with and 
without helical effects included, are performed against 
DNS computations. The spectral, statistical and spatial 
behaviors at large and intermediate scales of DNS flows 
are well restaured in both modeled flows. We notice some 
advantages for the model including helical effects, in par- 
ticular, concerning the evaluation of the helicity spectra 
and the probability distributions of the lateral velocity 
field gradients. LES PH also predicts a 4/3 spectrum 
for the helicity all the way to the cut-off, a point that 
will need further study. However, in our approach, we 
need to calculate at each time step the non local energy 
(and helicity) transfers, with an increased computational 
cost (roughly by a factor 2), when compared with other 
spectral LES models 0, . Whereas the role of helicity 
in fluids is not necessarily dynamically dominant, such is 
not the case in MHD where both kinetic and magnetic 
helicity play a prominent role, the former in the dynamo 
process and the latter in its undergoing an inverse cas- 
cade to large scales; furthermore, in MHD, spectra are 
not necessarily Kolmogorovian and the models presented 
in this paper may be of some use in MHD as well. The 
extension to MHD turbulent flows, with coupled velocity 
and magnetic fields, presents no major difficulties and is 
under study. 
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APPENDIX A: CLOSURE EXPRESSIONS OF 
TRANSFER TERMS 

For completeness, we recall here the expressions of the 
nonlinear transfer terms for the energy and the helic- 
ity, SE(k,p, q,t) and SH(k,p,q,t) respectively, under the 
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EDQNM closure assumption [l5[ : 

f E (k,t) = J J e km (t)S E (k,p,q,t)dpdq , (Al) 

T H (k,t) = jjj km {t)S H {k,p,q,t)dpdq. (A2) 

where A is the integration domain with p and q such that 
(k,p, q) form a triangle, and 6 k (t) is the relaxation time 
of the triple velocity correlations. As usual 9 h (t) is 
defined as : 

e kPq (t) = — — — , (A3) 

Pk + fJ-q + Mp 

where p,k expresses the rate at which the triple correla- 
tions evolve, i.e. under viscous dissipation and nonlinear 
shear. It can be written as: 

p k = uk 2 + \( y J^ q 2 E(q, t)dq) V2 . (A4) 

Here A is the only open parameter of the problem, taken 
equal to 0.36 to recover the Kolmogorov constant Ck = 
1.4. The expressions of S E (k,p,q,t) and SH(k,p,q,t) 
can be further explicited (with the time dependency of 
energy and helicity spectra omitted here) as : 

S E (k,p,q,t) = -b[k 2 E(q)E(p)-p 2 E(q)E(k)] 
pq 

- -£-c [k 2 H(q)H(p) - p 2 H(q)H(k)] 

= S El (k,p,q,t) + S E2 (k,p,q,t) 

+ S Ea (k,p,q,t) + S Ei (k,p,q,t) . (A5) 



Here, S El (k,p,q,t), S E2 (k,p,q,t), S E3 (k,p,q,t) and 
S Ei {k,p,q,t) are respectively used to denote the four 
terms of the extensive expression of S E (k,p,q,t). Note 
that S Es (k,p,q,t) and S Ei (k,p,q,t) are absent in the 
fully isotropic case (without helicity) , and, of course, all 
(k,p, q, t) terms below. 

S H (k,p,q,t) = -b[k 2 E( q )H(p)-p 2 E(q)H(k)] 
pq 

- -c[E(p)H(q) - H(q)E(k)} 
pq 

= S Hl (k,p,q,t) + S H2 (k,p,q,t) 

+ S H3 (k,p,q,t)+S Hi (k,p,q,t), (A6) 

with analogous short notations as before. 

In Eqs. (|A5[) and (|A6[) , the geometric coefficients 
b(k,p,q) and c{k,p,q) (in short, 6 and c) are defined as: 

b= P {xy + z z ): c= | z (i_ y 2) i (A7) 

where x, y, z are the cosines of the interior angles oppo- 
site to k, p, q. 

Let us now introduce a cut-off wavenumber k c , and de- 
fine the following three zones for the integration domain 
A of Eqs. (|Al|) and (|A2j) : the inner zone A< (with k, p, 
and q all smaller than k c , the cut-off wavenumber), the 
buffer zone A > (with p and/or q between k c and 3fc c ), 
and the outer zone A >> (with p and/or q larger than 
3fc c ); then, the boundaries of the transfer term integrals 
have to be adapted. 

In the inner zone corresponding to the fully resolved 
flow, the resolved transfers write: 



T£(k,t) = / -zP Q/3 (k)fc 7 / w ( g(p,t)u 7 (k-p,t)v a (-k,i)dpdk 

J\U\<k JO 
J|k|<fc JO 



(A8) 
(A9) 



r 



where P a p(k) = 6 a p — k a kp/k 2 is the projector on zone and the inner zone read: 
solenoidal vectors, as stated before. 

The transfers of energy and helicity between the buffer 
zone and the inner zone become: 



T>(fc,i) 
f>(k,t) 



3k c rk-\-p 



) kpq {t)S E (k,p,q,t)dpdq (A10) 



k—p 
3k c pk+p 



k—p 



kpq {t)S H (k,p,q,t)dpdq (All) 



and the transfers of energy and helicity between the outer 



T^(k,t) 



T>>(k,t) = 



oo pk-\-p 



3k(J k—p 
oo rk-\-p 

/ < 

3fc<y k—p 



kpq {t)S E {k,p,q,t)dpdq (A12) 
e kpq (t)S H {k,p, q ,t)dpdq (A13) 
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APPENDIX B: NUMERICAL 
IMPLEMENTATION OF THE MODEL 

As a first step, the Navier-Stokes equation is solved 
with the eddy viscosity and the helical eddy diffusivity 
v(k\k c ,t) and v{k\k c ,t) (see Eq. ([25]) . At this interme- 
diate stage, we obtain a partial estimation of the time 
updated velocity field, since the emission transfer term 
are not yet taken in to account. From this intermediate 
velocity, we compute the corresponding energy and he- 
licity density fields, say at wavevector k. They are then 
corrected with the appropriate EDQNM emission terms 
according to Eqs. (|2l))) and (|27p. In the time advance, 
the next step consists in reconstructing the three veloc- 
ity components from these updated energy and helicity, 
£(k, t) and H(k, t) respectively. When the velocity com- 
ponents are expressed as v a (k,t) = p a (k, t)e l ^ a ^ i ' t ' , the 
incompressibility condition leads to the following system 
of equations for the velocity phases and amplitudes (in 
short <j) a and p a ): 

P2P3COs((j>23) = 

pip 3 cos(4> 31 ) = 

PlP2COs((j)i2) = 
P2P3Sin((f>23) = 

Pip 3 sin(<j) 31 ) = 

■ U ^ fc 3^( k ) 

PiP2sm((j)i2) = — — 2 — (Bl) 



k 2 2£{k) 


-Pl{kl + kl) 


-p§(*?- 






2k 2 k 3 






fc|2£(k) 


- p\{k\ + k\) 


- P m- 






2kik 3 






fc§2£(k) 


-p\{kt + kl) 






2fcifc 2 



kfH(k) 

k 2 

fcltt(k) 

k 2 
k 2 H(k) 



with phase differences </> a /3(k, t) = 0^(k,t) — </> Q (k,i), a 
and P standing for the component indices. Note that 
when one component of the vector k is equal to zero, 
this case is treated separately in the code. Since only 



four of these equations are independent (because of the 
incompressibility condition), we are led to give an arbi- 
trary value for two of the variables. However, the choice 
of these arbitrary values is constrained. Indeed, from the 
set of equations Eqs. (|B1|) . we can derive an existence 
condition on the p a -amplitudes depending on the realis- 
ability condition (|7i(k)| < k£(k)) and which reads (with 
k-dependency omitted): 



k 2 



£(i-r) </£<(i 



7.2 
k 2 



£1 



(B2) 



with r = yl — H? jk 2 £ 2 . Thus, the p Q -amplitudes can 
be expressed as: 



pi 



Pa 



k 2 

W 1 



r» 2 



k 2 



(B3) 



where the i superscript denotes quantities based on the 
intermediate velocity field, which is a solution of the mod- 
ified Navier-Stokes equation Eq. (f2"5| . with eddy viscos- 
ity and helical eddy d iffusivity incorporated, and with 
r* = yl — TL i2 /k 2 £ l2 . Eq (|B3|) represents a projection 
of p l a (computed from the intermediate velocity field and 
which depends on £ % and Tt l ) to obtain the amplitude 
p a (k, t) (depending now on £ and H.) at the updated 
time step. This allows not to modify the velocity field 
when f % g = and = 0. 

If one or two componants of the k- wavevector are equal 
to zero, the set of equations (|B1[) is rewritten from the 
divergence-free condition. Apart from this, the recon- 
struction procedure is similar. 

Finally, to rebuild the different velocity phases, <j)\ is 
assumed to be fixed to its value given by the intermediate 
component v\ (k, i) . The set of equations (|B1[) is then 
solved and we obtain the updated Fourier velocity field. 
Note that a different choice for the fixed phase leads to 
no significant changes in our numerical tests. 
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